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Abstract. We consider the semiclassical limit for the nonlinear Schrodinger 
equation. We introduce a phase/amplitude representation given by a system 
similar to the hydrodynamical formulation, whose novelty consists in including 
some asymptotically vanishing viscosity. We prove that the system is always 
locally well-posed in a class of Sobolev spaces, and globally well-posed for a fixed 
positive Planck constant in the one-dimensional case. We propose a second order 
numerical scheme which is asymptotic preserving. Before singularities appear in 
the limiting Euler equation, we recover the quadratic physical observables as well 
as the wave function with mesh size and time step independent of the Planck 
constant. This approach is also well suited to the linear Schrodinger equation. 



1. Introduction and main results 

1.1. Motivation. We consider the cubic, defocusing, nonlinear Sclirodinger equa- 
tion (NLS) 

ie^tu^ yAit'^ = |n^pu^ {t,x) £R+ xR'^, (1.1) 
with WKB type initial data 

n"(0,x) = ao(x)e^<^o(^)/", (1.2) 

the functions oq and (j)Q being real-valued. The aim of this article is to construct an 
Asymptotic Preserving (AP) numerical scheme for this equation in the semiclassical 
limit e — 7- 0. We seek a scheme which provides an approximation of the solution 
with an accuracy, for fixed numerical parameters At, Ax, that is not degraded 
as the scaled Planck constant e goes to zero. In other terms, such a scheme is 



consistent with (1.1) for all fixed e > 0, and when e — )• converges to a consistent 



approximation of the limit equation, which is the compressible, isentropic Euler 



equation (1.4). The difficulty here is that, when e is small, the solution u'^ becomes 
highly oscillatory with respect to the time and space variables, and converges to its 
limit only in a weak sense. In order to follow these oscillations without reducing At 
and Ax at the size of e, which may be computationally demanding, our construction 



relies on a fluid reformulation of (1.1), well adapted to the semiclassical limit. 

In kinetic theory, plasma physics and radiative transfer, the interest in AP 
schemes has considerably grown in the recent years since the pioneer works on 
the subject |19| l26l \3T\ 136) . Among a long list of works on this topic, we can men- 
tion [Ml |28l [III IMl El [2Ql [251 [IS 1211 |38l I23]- For the Schrodinger equation, fewer 
works exist. In the stationary case, one can cite |3] which proposes a WKB-type 
transformation in order to filter out the oscillations in space. In the time-dependent 
case, the usual numerical methods for solving the nonlinear Schrodinger equation - 
finite-difference schemes |5H [2^ [H [3l] , time-splitting methods |46| \50\ [9] , relaxation 
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scheme [8j or even symplectic methods pSj - are designed in order to guarantee the 
convergence of the discrete wave function when e > is fixed. As it is analyzed 
for instance in |44| H5] for the case of finite-differences or in [37] for Runge-Kutta 
schemes, these methods suffer from the oscillations in the semi classical regime if the 
time and space steps are not shrunk. Note however that, among these methods, the 
status of splitting methods is particular. Indeed, in the linear case, it is shown in 
|5] (see also the recent review article [32]) by Wigner techniques that, if the space 
step Ax has to follow the parameter e when it becomes small, the time step can be 
chosen independently of e (before the appearance of caustics) if we are only inter- 
ested in getting a good approximation of the quadratic observables associated to . 
Still, none of these methods is Asymptotic Preserving in the sense defined above. 

The quantum fluid reformulation of the Schrodinger equation provides a natural 
framework for the construction of AP numerical methods. The Madelung transform 
|42) is a polar decomposition of the solution of (1.1), written as 

^^(t,^) = a^(^,x)e^*'(*'^)/^ 

where the amplitude a"^ and the phase 0^ are real-valued. Inserting this ansatz into 
(1.1 1 yields the system of quantum hydrodynamics (QHD), or Madelung equations, 
satisfied by = (a^)^ and = V0^: 

' dtp' + div ip'v') = 0, p^^^g = (ao)2, 

dtv + ■ Vv' + Vp' = -V [-^) . y\t=o = V^/'o. 

This system takes the form of the compressible Euler equation (with the pressure 
law p{p) = p2/2), with an additional term testifying for the quantum character 

of the equation, the so-called Bohm potential ( ^^f- )■ Two comments are in 
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order. First, the term in e does not appear any more as a singular perturbation 
and, in the limit e — )• 0, the quantum pressure disappears, leading formally to the 
compressible, isentropic Euler equation 

(dtp + diY{pv) = 0, pi^^o = (ao)2, 

which is the semiclassical limit of our problem (actually valid before the appearance 
of shocks). Second, the fluid-like form of this system enables to consider many 
numerical methods originated from computational fluid mechanics. 

For the linear Schrodinger equation, the Madelung formulation is at the heart 
of the method of quantum trajectories (or Bohmian dynamics), see |52j . which are 
based on a Lagrangian resolution of this system. In Eulerian coordinates, the work 
|22] also exploits this formulation to construct Asymptotic-Preserving schemes for 
(1.3). Unfortunately, the main drawback of the QHD system is the form of the 
quantum potential, which becomes singular as the density p^ vanishes (see [13] for 
a recent survey). Hence, these methods do not provide AP schemes in the presence 
of vacuum. Another point of view was recently developed in |15) for the nonlinear 



equation (1.1 1, taking advantage of another fluid reformulation of the equation - due 



to Grenier |29j - which tolerates the presence of vacuum. In the next Subsection 1.2 



we present in detail this formulation relying on a complex amplitude version of 
Madelung transform. However, as we will see, the drawback of this reformulation is 
that this system may develop singularities even for fixed e > 0, whereas the solution 
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of NLS remains smooth. Then, in the last Subsection 1.3 of this introduction, we 
present a modification in order to remedy this problem and preserve the regularity 
of the solution. The new QHD model that we propose contains a viscosity term in 



the equation for the velocity, but is still equivalent to the original equation (1.1). 



1.2. Backgrounds on NLS and its semiclassical limit. In this section, we recall 



some known results on the Cauchy problem for NLS (1.1) and on its semiclassical 
limit as e — 7- 0. For details on this subject, we refer to the textbooks |17| and |12) . 

The following result is standard in the case e = 1, and can easily be deduced 
when e G (0, 1], by scaling arguments. 

Proposition 1.1. Letd(^3, and G H'^{R'^). 

(i) There exists a unique solution G C(M; H^{W^))nL\l^{M.; ^^^'^(M'^)) to ( [LTj ), 
such that tijj^Q = Uq. Moreover, the following conservations hold: 

Mass: ^^Wu" {t)\\\2 = ^ . 
d f 

Momentum: "i- Im / u^{t, x)'S/u'^ {t, x)dx = 0. 
dt J^d 

Energy: ± {\\eVu%t)\\l, + \\u%t)\\l) = 0. 

(ii) If m addition u% G H^{W^) fork(^n,k~^ 2, then G C(M; H^{R'^)). 
Note that, if = a^e**^"/^ with (f)^ G M, the above three conservation laws become 



±\\a%t)\\l,=0. (1.5) 



^ ( Im [ a''{t,x)\/a%t,x)dx +- [ \a^{t,x)\'^\/(j)%t,x)dx] =0. 

I (||£Va^(t) + za%t)V^%t)\\l, + ||a^(t)||i4) = 0. (1.6) 



Let us now describe the limit e — )• for (1.1). The toolbox for studying semi- 
classical Schrodinger equations contains a variety of methods, depending on the 
quantities we are interested in. If one is only interested in a description of the 
dynamics of quadratic observables, such as the mass, current or energy densities, 
the Wigner transform is well adapted and has been applied to linear Schrodinger 
equations and Schrodinger- Poisson systems \27 \ I4 H |55]: on the other hand, it is not 
adapted to the study of NLS |14) . Still for a description of quadratic observables. 



the modulated energy method \10\ |5^ HOl [2] enables to treat the nonlinear equation 



(1.1 ). Yet, for a pointwise description of the wavefunction u"^, WKB techniques are 
more convenient. We refer to |12| for a presentation of WKB analysis of Schrodinger 
equations. 



In the case of (1.1), also referred to as supercritical geometric optics in this 
context, the justification of the WKB expansion was done by Grenier |29) . Let us 
briefly present his idea. As we said above, the main inconvenience of the QHD 



system (1.3) is the singularity of the Bohm potential at the zeroes of the function 
{i.e. at the zeroes of the wavefunction u'^). Looking again for solutions under 
the form 

7x^(t,x) = a^(t,x)e^'^'(*'^)/", (1.7) 
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but allowing now to take complex values, Grenier proposed to define and (f)'^ 
as the solutions of the system 



2 



0, 



*|t=o 



'|t=0 



[1.8) 



This system, which is formally equivalent to (1.1), can also be expressed in terms 
of the unknowns and = Vcp^. Differentiating with respect to x the second 
equation in (1.8), we find, using the fact that is irrotational: 



1 £ 

dta" + • Va*" + -a^divz;^ = i-Aa^, 



e\2 



0, 



:i.9) 



The important remark made in |29j is to notice that the above system is hyperbolic 
symmetric, perturbed by a skew-symmetric term. In the limit case e = 0, and if oq is 
real- valued and nonnegative, (1.9) is the Euler equation (1.4) written in symmetric 
form (see e.g. j43l 118) ): 



dta + V ■ Va -|- -a div u = 0, 
dtv + vyv + V\a\'^ = 0, 



a\t=o 

V\t=Q 



Oo, 

V0O. 



:i.lO) 



This remark enables to obtain, in small time (but independent of e G [0, 1]), uniform 
estimates of the solution {a^ ,v^) of ( |1.9[ ) in Sobolev norms and prove the convergence 
of these quantities to the solution {a,v) of (1.10), as long as this solution exists and 
is smooth. From the numerical point of view, the regular structure of (1.9), in 
particular when vanishes or has very small values, is an advantage compared to 
(1.3), and is at the basis of the schemes constructed in |15j . 



1.3. New model and main results. A drawback of (1.9) is that large time ex- 
istence results (for fixed e > 0) do not seem available. In particular, it is not 



known whether (1.9) still has a smooth solution past the time where the solution 
to (1.10) has ceased to be smooth (such a time necessarily exists if and oq are 
compactly supported, regardless of their size, from |43| — see also [53j). A practical 



consequence of this fact is that a numerical method based on (1.9), such as the 
one proposed in |15| , may not approximate correctly the original equation (1.1) on 
arbitrary time intervals, for fixed e > 0. 



To overcome this issue, we use the degree of freedom given by (1.7) in a manner 



which is slightly different from the approach introduced by E. Grenier. From (1.7) 
we have 



1, 



q2 \ rfiJ4>^le 
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Introduce the viscous term aF e^'^^ ^ ^ Acj)^ , and reorder the terms so as to consider 
the new system 



dtcP' + ^\Vcl)'\^ + \a'\^ = e^Acl)', 
Proceeding hke above, we get the following new system: 

dta^ + v"^ ■ Va^ + -a^ divv^ = i-Aa^ — iea^ divf^, 
2 2 

dtv' + -Vv" + 2 Re (a^Va^) = e'^Av", 



:i.i2) 



We recover physical observables, respectively particle, current and energy densities, 
thanks to the following formula 



P 

y = elm(a^Va^) + p'^v'^, 
Remark 1.2 (Linear case). In the case of a linear Schrodinger equation 



:i.i3) 



-Au' = Vu" 



we will see that the introduction of this artificial diffusion is even more striking than 



in the nonlinear framework; see Subsection 4.1 



Remark 1.3 (Irreversibility). Our choice for introducing the viscous term in ( |1.11[ ) 
makes the system irreversible (in time), while (1.1) is reversible. To solve the 
Schrodinger equation backward in time, it suffices to change the signs to 



-e^A^^ 



'\t=o 



Let us now present our main theoretical results on this system. Our first result 
concerns the local in time well-posedness of the Cauchy problem in any dimension, 
for fixed values of e ^ 0. 

Theorem 1.4. Let (ao,(Ao) e H%R'^) x H'+'^iR'^), with s > d/2 + 1. Then the 
following holds. 

(i) The system (1.12) admits a unique maximal solution [a'^ ,v'^) £ C{[0,T^^^); H'^ x 

(ii) There exists T > independent of £ £ [0, 1] such that T^^x ^ ^- Moreover, the 
L°°([0, T]; ff'' X H'^) norm of {a^,v^) is bounded uniformly in e € [0,1]. 

(ill) Defining 0^ G C([0, T^^J; i?«+i) by 



(f>^t,x) = (f>o{x) 
then the function 



\v'^{s,x)\'^ + |a^(s,2;)p - div v'^{s,x) \ ds, 



:i.i4) 



is the unique solution to (1.1) satisfying (1.2). 
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Remark 1.5. In fact, the proof of this theorem provides a criteria of global exis- 



tence for the solution to (1.12). Indeed, we shall prove that the lifespan T^^x 



independent of s > d/2 + 1 and that 



^ max 



(a^7;^)(^)||vKl.- + \\a^t)\\l^) dt 



oo. 



:i.i5) 



Our second result states the convergence of the solution to (1.12) towards the 
solution of the Euler equation, as e — )• 0. 

Notation. For two positive numbers and the notation < means that 
there exists C > independent of e such that for all e G (0; 1]) ck"^ ^ 

Theorem 1.6. Let s > d/2 + 3, and (ac^o) € H'^ x H'+^ . LetT > such th at the 
Euler equation (1.4) has a unique solution G C{[0,T]; H'^ x H^). Then (1.10) 

has a unique solution {a,v) G C([0, T];if** x H^). Moreover, for e > sufficiently 
small, we have Tii ^ T and 



\\{a^,v'') - (a,';;)||ioo([o,T];H- 



Theorems 1.4 and 1.6 are not different from the corresponding results obtained by 



Grenier in [29] on his system (1.9). We now state a global existence result for fixed 



e > in dimension one, that is not available on (1.9). This result takes advantage 



of the viscous term in the second equation of (1.12) 



Theorem 1.7. Suppose d 



1 and let (ao, <po) G x H^'^^ , with s ^ 2. Then for 

oo in 



e > fixed, the solution to (1.12) is global in time, in the sense that 
Theorem 11.41 



'£ 

max 



This article is organized as follows. Section [2] is devoted to the proofs of our 
three theorems: the local existence result in Subsection 2.1 the semiclassical limit 
in Subsection 2.2 and the global existence result in Subsection |2.3| Section [3] is 
of different nature and concerns numerics. We describe a second-order Asymptotic 
Preserving numerical scheme for NLS, based on the reformulation (1.12), and we 
give the results of numerical experiments, in dimensions 1 and 2. We show that 
our scheme is indeed AP before the formation of singularities in the Euler equation, 
which means that the time and space steps At and Ax can be taken independently 
of e. After the formation of singularities, we show that, in order to get a good 
approximation of the solution, we need to take At and Ax of order 0(e). At the 
end of the article, in Section]!] we sketch two extensions of our results: we discuss on 
the case of the linear Schrodinger equation, then we discuss on other nonlinearities. 



2. Proofs of the main theorems 



This section is devoted to the proofs of our main results stated in Theorems 1.4 
[LHlandOl 



2.1. Local existence. In this paragraph, we show the local existence of a smooth 
maximal solution to (1.12). 
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Proof of Theorem 1 A , Items (i) and (ii): local existence and uniform estimates. 
Following |29) . introduce the vector- valued unknown function 

/Rea"\ /af\ 



Im 



pd+2 



In terms of this unknown function, (1.12) reads 



dtU' + ^A^{U')djW = ^LU' + e^DU' +e^B^{U')djU', (2.1) 

where the (fc, ^)isgfc,^sgd+2 elements of the matrices A^{U) G are given by 

f K,ki^) = Uj+2 for k = l,...,d + 2 

Al^^,iU) = Ui/2, A{^^,iU) = U2/2 

A^:+2,i(f/) = 2t/i, A^^,^,{U) = 2U2 

Aj, ^(U) = otherwise. 

The linear operator L is given by 

/O -AO ... 0^ 
L = A ... 

\0 Odxd 

A first important remark is that even though L is a differential operator of order two, 
it causes no loss of regularity in the energy estimates, since it is skew-symmetric. 



S 



Also, (2.1) is hyperbolic symmetric (or symmetrizable). Indeed, let 

72 

This matrix is symmetric and positive, and SA^ is symmetric for all k 

SA^{u) € cSd+2(M), yu e 

The diffusive term D is given by 

/O ... 0^ 

L> = ... 

\0 IdxdA 

Finally, the matrices , defined analogously to A^ , are given by 
r Bl^^,{U) = U2, Bl^^,{U) = -Ui 
[ B-j^ ^{U) = otherwise. 



(2.2) 



When the right-hand side of (2.1 ) is zero, local existence of a unique solution in H'^ 



with s > (i/2 -|- 1 is a consequence of standard quasilinear analysis for hyperbolic 
symmetric systems; see e.g. j49j . The important point in energy estimates consists 
in computing 

i- {A'{SU'),A'U') , A = (1 - A)^/\ 
dt 
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and taking advantage of symmetry to obtain, thanks to tame estimates, 
d 



dt 



for some constant C depending only on A and d. 

In the case of the complete system (2.1), we first recall that as noticed in 
the term L does not alter the above energy estimate, since SL is skew-symmetric, 
so we have 

^{A'{SU'),A'U') ^ CWU'Wwi,'^ {A'{SU'),A'U') +2e^ {SA'DU',A'U') 

d 

+ 2e^{SA' {B^{W)djU') ,A'U') . 
i=i 

In terms of (v^ , a"^), the first new term reads 



2e' {SAWU^,AV^) = y E/ ^i^'vDA^vl = / 1 

k=l 



VA'v 



For the second new term, we have 
2e{SA' [B\U')djW) , A'U') = 



A' {a'^djv^j) A'al 



A' (aldjv'j) A'al 



By Kato-Ponce commutator estimates |35) . 

||A" {aldjv'j) -alA'djv'jW^^ < || Va|||L- llt-IH//- + ||a|lkH|Vt;^IU- • 
Then using Cauchy-Schwarz and Young inequalities. 



alA'djv'A'al 



Choosing 5 > sufficiently small and independent of e, we infer 



d 



£ 



{A'{SUn,A'U') + ^||VA^^;^||i. < ||C/' 



iH" "I" iI'^^IIl°° 11'^'' ■ 

(2.3) 

The local existence result then follows easily by adapting the standard arguments 
from |49) . In view of the energy estimate (2.3), we also infer (1.15). The lifespan 
^max of the maximal solution to (1.12) must be expected to actually depend on 
e. For instance, if (t)o,ao G Cg°(M''), then from |43j , T^^x < regardless of the 
dimension d. On the other hand, we prove further (proof of Theorem 1.7) that if 
d 



1 

) max 



oo for all e > 0. 



Item (in): is solution of the original NLS problem. Define 0^ by (1.14). Let 
us prove that (a^,0^) G C{[0,T^^^); x solves the system (1.11). Since 

(a^ v') G C([0, T^^^); H'xH'), we readily have cp' G C([0, T^^J; H'-^). Moreover, 
is irrotational: indeed, curlu^ satisfies a homogeneous equation (see e.g. jJl 
p. 291]), it is zero at time t = 0, so curlv^ = 0. Set 17^ = Dv^ — Vu^, where Dv^ 
stands for the Jacobian matrix of v^, and Vv^ stands for its transposed matrix. It 
solves 

As a consequence, V|f^p 



second equation of (1.12) 



2v^ ■ Vv^. We then deduce from (1.14) and from the 

dt -v') = vdt^' - dtv' = 0, 
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SO we infer from the initial condition v^^^q = Vcpo = "^(p^^^Q that v'^ = V^^. There- 



fore, V(A-EC([0,r^,J;F^ 



Replacing with V(/)^ in (1.14) shows that (j)'^ solves the first equation in ( |1.11[ ), 
and the second equation in (1.11) is obviously satisfied. 



Now, define by (1.7). The property u'^ S C([0, T^]; i?**) is straightforward. 



since H'^iW^) is an algebra. It is clear that solves (1.1) and (1.2). Since 
C([0, T'];H\M!^)) C C([0, T']-H\M.'^)) n lS/'^([0, T"]; VF^'^CM"')), 
uniqueness stems from Proposition |1.1| □ 

2.2. Convergence towards the Euler equation. In this subsection, we study 



the semiclassical limit of (1.12). 



Proo f of Theorem 1.6, The proof proceeds as for [29| Theorem 1.2]. First, The- 
yields a solution (a, f ) G C([0, Tj^^x); -f^* ^ 



orem 



1.4 



Necessarily, T^ax > ^ 



(where T is an existence time for the Euler equation) , for if we had T^ax ^ ^) then by 
uniqueness for (1.4), (|ap,t;) = (p, f) G L°°([0, T]; iJ* x i/*). The equation for a in 
dTlUl ) then shows that a £ L°°([0, T]; i?*"!), hence {a,v) G L°°([0, T]; if'-^ x H'), 
which is impossible if T^^^x ^ ^- this yields the first part of the theorem. 

To prove the error estimate, introduce the vector-valued unknown U, associated 
to ia,v). We know that U G C{[0,T]; H'). Consider the error = - U. It 
solves 



Rewrite this equation as: 



3=1 



The operator on the left hand side is the same operator as in (2.1 ). The term LW 



is not present in the energy estimates, since it is skew-symmetric. The term eLU is 
considered as a source term: it is of order e, uniformly in L°°([0, T]; i/**"^). Next, 
the first term on the right hand side is a semi-linear perturbation: 



\{A^{U^)-A^{U))d,U\\^,^, ^ \\{Aim-A^{U))\\ 



H 



s-2 



\u\ 



<\\[A^{W' + U)-A\U))\\^, 



s-2 

L^,\\U Wh"-'^) II ~ ll_H"''-2) 



where we have used tame estimates. Finally, we know that is bounded in 
L°°([0, T] X M'^), as the difference of two bounded terms. The terms involving D 
and B are treated in a similar way, by adapting the estimates used in the proof of 
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Theorem 1.4 We end up with 



(ft 



<e^ + {Sk'-^W\k'~^W'), 
and the result follows from Gronwall lemma. □ 

2.3. Global existence result in one dimension. In this section, we consider 
(1.12) with e > fixed, in dimension d = 1. To prove a global existence result, we 
shall need finer estimates than (2.3), in particular for low order derivatives. 



Proof of Theorem 1.7 . For s = and s = 1, (1.12) can be refined, and this will be 
useful to prove a global existence result. 



estimates. Multiplying the second equation in (1.12) by a^, integrating in space, 
and considering the real part yields: 



du e 

dt" 



L2 



0. 



We have simply recovered the conservation of mass (1.5). 

Multiplying now the first equation in (1.12) by and integrating in space, we 
take advantage of the dimensional one: 

- — II 



\a^{t,x)\'^\dxV^{t,x)\dx. 



(2.4) 



estimates. Differentiating (1.12 ) in space, multiplying then the first equation by 
dxV^, the second by dxO,'^, and integrating, we find (thanks to the symmetrizer): 

d 



^ {\\dxv^mh+mam)\\i,)+e^dyit)\\i. 



ioo + \\a 



) 



m + II" 



2 ^ 



(2.5) 



Now, we have the tools to prove Theorem 1.7 In view of Remark 1.5 it suffices 
to prove that for any fixed e > 0, 

v',dxv',a',dxa',\a'\^ G ( [0, oo); L~(M'^)) . 

Since e > is fixed, we shall omit it in the notations, and fixed it equal to 1 in 
(1.12). Using Item (Hi) of Theorem 1.4 and in view of (1.6), we have the a priori 
estimate 

l|a(t)||L4 < 1, Vt^O. 
The refined estimate (|2.4|) for v and Cauchy-Schwarz inequality yield 



1 d 



2 dt 

Using Young's inequality 



vmh + \\dMmh 



< 



L4 



\dxv{t)\\L2 < \\dxv 



L2- 



\dMmL'^^s\\dMmh + 



1 

4^' 



we infer, for 5 > sufficiently small. 



d 

~dV 



h + \\dxvmi2<i. 
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Therefore, 



vmh+ f\\dMr)\\hdT<l + t. 
Jo 



10 

The Gaghardo-Nirenberg inequahty yields 

Jo Jo 

t 



<(l + t)i/4 / ^i + \\dMr)\\h)dT<{l + tf\ 





where we have used the Young inequahty. In view of ( [1.6^ , we infer 

\\dMt)\\L^ < 1 + HtHt)\\L2 < 1 + \\v{t)\\L^, 
and the Gaghardo-Nirenberg inequahty yields 

\Ht)\\L^ < \\a{t)\\];,'\\dMt)\\]!' <^ + \\vm]!i. 

We infer 

[\\Hr)\\L- + \\a{r)\\l^)dT<{l + tf/\ 
Jo 



It only remains to prove that d^v, d^a € ^[^^([O, oo); L°°(R)). The estimate ( [2^ 
yields 



^0 

t 



vir)\\hdT < 1 + / {\\dMr)\\L- + 1 + \Hr)\\l^) (||a(r)||?,i + Mt)^^^) dr 
Jo 

<l+ f {\\dMr)\\L- + 1 + \Hr)\\l^) (1 + \Hr)\\l^ + \Ht)\\\.) dr 
Jo 



< 1 + / [\\dMT)\\L^ + 1 + lla(r)lli^) (1 + \\v{t)\\1,) dr, 





from Sobolev embedding. Using an integration by parts and the Young inequality, 
we find 

\\dMr)\\h^\Hr)U4dlv{T)U. 

^ lia f Ml 1 1 I II ^ MI2 \\9xV{t)\\12 

\\dxV[T)\\L^ +l + \\a{T)\\l^ 



\\dMT)\\L^±J^+ l|a(T)||ioo„^,,^,„2 

id 



ll-^xi^V ^ll-'.'" ll^-VVIIL"" II / \ 



Taking 6 > sufficiently small, we infer: 

Mr)\\hdT <l + {l + t) f {\\dMr)\\L^ + 1 + lla(r)llioc) dr 



t rt 
\\dh 

JO 

rt 



+ {l + t) [ {\\dMr)\\L^ + i + MT)\\l^fdT 
Jo 

<l + t^ + {l + t) [ \\dMT)\\l^dT + {l + t) [ \\a{T)\\UdT. 
Jo Jo 
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We also have 



t ft rt 

4 J, ^ + , / ll„,^,M|2 



a{r)\\lo.dT<t+ \\v{T)\\i^dT<t+ \\v{T)\\L2\\dMr)\\L^dT 

JQ Jo 

<t+VlTt[ \\dMr)\\L^dT<{l + tf/'', 
Jo 

and so 

r \\dlv{T)\\l,dT < (1 + tf/' + {l + t) f \\d.Mr)\\lo^dT. 
Jo Jo 

Now Gagliardo-Nirenberg and Cauchy-Schwarz inequalities yield 

* \\dMr)\\lo.dT < f \\dMr)\\LA\dlv{r)\\L2dT 
Jo 

1/2 ^ rt \ 1/2 



< 



\dMr)\\hdT] ( / \\d-iv{T)\\i,dT 







< ((1 + tf'^ + + \\d.v{r)\\lo^dT 

<(l + tf/4 + (l + t) (^j\\dMr)\\l^dT 



1/2 

ir J 

1/2 



The Young inequality implies 

f\\dMr)\\l^dT<{l+t)\ 
Jo 

hence d^v G L2,,([0, oo); L^) C L1„,([0, oo); L^). 

Since O^^a G LJq^([0, oo); L^), Gagliardo-Nirenberg shows that it suffices to prove 
that d'^a G L\^^{\Q^oo); L'^) to conclude that G LJq^([0, oo); -L°°). For that, we 



use again Item ( in) of Theorem 1.4 from which we now that ae^* G L\^^{[Q,oo); H"^). 



Differentiating ae*'*' twice, we have, since dx(t> = v, 

\\9la{t)\\L2 < \\ae'''^\\H2 + \\vdxa\\L2 + ||a5a;v||i2 + ||u^a||L2 

< ||ae*''^||j|^2 + ||t'||Loo||9a,a||i2 + ||a||Loo ||52,t;||i2 + ||w||ioo. 

From the above estimate, we infer d'^a G -^^ioc([0, oo); L^), which completes the proof 
of the theorem. □ 

Remark 2.1 (Higher dimension). Even though we have obviously taken advantage 
of the one-di mens ional setting to use the embedding C L°^, the true reason 
why Theorem 1.7 is limited to d = 1 lies elsewhere. The refined estimate (2.4) 



becomes, if ^ 2, 

1 d 
2dt"^ 



l2+e^Vv^t)\\l2< I {\a^it,x)\' + \v^it,x)\^)\Vv%t,x)\dx. 



Note the appearance of the new term J \v^{t, a;)p|Vf^(t, x)\dx: it correspond to the 
fact that J ■ {v^ ■ Vv^) is zero if d = 1, but not if d ^ 2 (in general). This aspect 



seems to be the only real limitation to extend Theorem 1.7 to d ^ 2. 
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3. Numerics 



In this section, we define a second order numerical scheme for ( 1.12 ). Since e is no 



longer a singular perturbation parameter in this reformulation of NLS, this scheme 
is naturally Asymptotic Preserving for our original problem in the semiclassical 
limit e — 7- 0, as long as the solution of the Euler equation remains smooth. Since 
we solve our problem in dimension 1 on a bounded interval and in dimension 2 



on a rectangle, we add periodic boundary conditions to (1.1) or to (1.12), both 
formulations remaining equivalent. 



3.1. The AP numerical scheme. The nonlinear system to be solved reads 



dta + div(at;) + {ie — -)a div v = «-Aa ; O|t=o = '^Oi 

2 

dtv + V + \a\^j = e^/^v ] v\t=Q = '^o- 



(3.1) 



The semi-discretization in time is realized through a Strang splitting scheme. Let us 
denote by At„ = t„ — t„_i the variable time step, such that t„ = Yllt=i ^^k- On each 
time step [t", t""*"^], we split (3.1) into two subsystems and apply the second-order 
Strang splitting algorithm. 



Step 1 for t G [tn, tn + 




2 



h\t=t„ 



Step 2 for t G t„ -|- At„+i]: 



dta2 + div(a2W2) + {ie - -)«2 divw2 = ; a2\ 



dtV2 + 



-n — ^ "2 



; 



t — tfi 



V2\t=t„ 



n + l ) 
2 



V , At 
l\t=t„ ■ 



n + l • 



Step 3 for t G t„ -|- 



dm = e'^Av3 ; V3it=t„ = V2\t=t„+^- 



Then, (03,^3) _ At„+i is an approximation of {a,v)u^i , solution to (3.1). 

1^ — in~\ 2 ^ 



This splitting scheme enables to decouple the fluid part (step 2) from the parabolic 
and Schrodinger parts (step 1 and 3), which allows a standard treatment of each 
step. Note that for e = 0, we recover exactly the Euler equations ( |1.10 ). We denote 
by (a"'^, w""'^) the numerical solution at time i„. In the spirit of |30) . the following 
result can be proven. 



Proposition 3.1. Under the assumptions of Theorem 1.4, there exist Eq > and 
C,cq independent of £ ^ [0,eo] such that for all n & N such that tn G [0,r], where 
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T is as in Theorem 1.4 (ii), for all At„ € (0, cq], 



a 



eIm(a"'^Va"'^) + \a 



Li(Rd)nL°°(Rd) 



n,£\2 yn,e 



^ C(maxAt„)'', 



^ C(maxAt„) 



We assume, without loss of generality, that we are in a periodic framework for 
the space variable. This allows to solve steps 1 and 3 in a spectral way by using fast 
Fourier transform, whereas for step 2, we use a Lax-Wendroff scheme (with direc- 
tional splitting in dimension 2). The extension to Dirichlet or Neumann boundary 
conditions would require to consider sine or cosine transforms in place of FFT. The 
total scheme is consistent with ( |3.1[ ), and second-order in time and space. Since the 
Lax-Wendroff scheme is explicit, a Courant-Friedrich-Lewy condition is necessary 
for stability: 

Ax 



Atr, € CFL- 



max,- 



+ \a]\ 



where denotes the approximation to v{t„ 



being a node of the mesh 



discretization of the computational domain Q, and Ax is the space discretization 
parameter. In all our numerical experiments, we take CFL = 0.8. 

3.2. Numerical experiments in dimension 1. In order to analyze the behavior 
of the numerical solutions (3.1), we proceed like in |i6j, and compare the behavior 



of numerical physical observables (1.13) computed by our AP numerical scheme 
to a reference solution. Since exact solutions to (1.1) with initial data (1.2) are 



not available, we numerically generated one thanks to the usual splitting method 



described in [6], applied to (1.1 )-(1.2 ) with a meshing strategy that ensures that the 



time step At < e and the space step Ax < e. We consider the initial datum 

ao(x) = e-25(-o.5)^, 

vo{x) = -^d^ In (e5(^-0-5) + e-5(x-o.5))^ 

with X G [—1/2, 3/2]. In all our simulations, the smallest value of e is 10~^. We then 
discretize the x-interval with J = 2^^ subintervals and the time step is At = e/100. 
We refer to the reference solution in a generic way by the notation n^^j. With the 
considered initial datum, new oscillations appear in the reference solution passed a 
time r* ~ 0.10, meaning that singularities were created in the limit Euler system 
(1.10). In Figures [T] and |2j we plot particle and current densities p"^^^ and j^^j for 



s = 10 respectively at time t = 0.05 and t = 0.13, before and past formation of 
shocks. We present on Figure |3] a zoom close to singularities which shows oscillations 
on physical observables. 

In order to show the asymptotic preserving property of our scheme, we compute 
relative £^ errors at time t„ thanks to the formula 



err„£ 



(tn) 



errje (t^ 



llPref 
lUrcf " 



P 



1, 



• n,e 



I n,£\2 

a- Va 



'3 > 



where the £^ norm is 



mi 



Ax 



j-i 

k=l 
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(a) Particle density prof 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



(b) Current density jrof 



Figure 1: Particle and current densities at time T* = 0.05 and for e = 10 ^. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



(a) Particle density prc^ 



>-0.05 
'-0.1 

-0.15 




0.2 0.3 0.4 0.5 



(b) Current density jrcf 



Figure 2: Particle and current densities at time T* = 0.13 and for e = 10 ^. 



We evaluate the error function for different values of J = 2*^, where M is an 
integer here chosen in [4, 13], and various scaled Planck parameter e in [10-^10°]. 
We plot on Figures |4) [5] [6] and [7] the relat ive error on physical observables p'^ and 

computed with our AP schemes at time t = 0.05 and t = 0.13 respectively before 
and after the formation of singularities. We clearly see that the error is proportional 
to (Ax)^, independently of e before formation of singularities. Indeed, the mean 
slope of the error curves with respect to A.x is close to the value 2 (subfigures (a)). 
The independence with respect to e appears in subfigures (b) where error curves 
are flat. After the formation of shocks, the behavior of our scheme is very different 
and the mesh parameters have to be reduced as e — )• to get good accuracy. 

For a comparison, we present the same study for the classical time splitting 
scheme applied to (1.1 )-(1.2 ) on Figures [8) |9j [To] [n] Thanks to the spectral accuracy 
of Fast Fourier Transform, the error levels are smaller than for our AP scheme which 
is only second order with respect to time and space variables. However, contrary to 
our AP scheme, the error curves clearly depend on e, the error being of order 0(1) 
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(a) Particle density 




0.63 0.64 0.65 0. 



(b) Current density jrcf 



Figure 3: Zoom around singularity of density and current at time T* = 0.13 and 
for e = 10-^. 




-e^e=1 .0000 
-4 -E=0.6400 
->- E=0.3200 
□ E=0.1600 

E=0.0800 

-* -E=0.0400 
-t^- E=0.0200 
o E-0.0100 
E=0.0050 
- :- -E=0.0025 
E=0.0010 
E-0.0001 



- ^- -<3- --^ - < - -c\~ ^ 



J=1024 
J=2048 



(a) Error w.r.t J = 2^* 



(b) Error w.r.t e 



Figure 4: errpe(t = 0.05) for AP scheme 



when e is smaller than 2.10 In order to have good accuracy, it is required to take 
Ax<e. 

To make our presentation complete, we present the reconstruction of thanks to 



equation (1.14) with e = 0.005. We therefore compute the phase cf)'' with amplitude 



a"^ and velocity u^. Since we have access to values of these quantities at every 
discrete time tn, we approximate the time dependent integral thank to a simple 
rectangular quadrature. We see that we can recover a very good approximation of 
the wave function with very few points before the formation of singularities (see 
Fig. 12). Obviously, one needs more points to have a good reconstruction of the 



wave function passed the singularities (see Figs. 13 and 14) 



AP SCHEME FOR NLS 



17 
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- E=0.0025 
-E=0.0010 

6=0.0001 



a: 10-", 





J=16 
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-\>- J=64 


□ J=128 




J=256 




- ft -J=512 




-+- J=1024 


^ at] 


o J=2048 




— 9— J=4096 


□ ^"''^ 


- e- -J=8192 




(a) Error w.r.t J = 2^^ 



(b) Error w.r.t e 



Figure 5: errje(t = 0.05) for AP scheme 
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(a) Error w.r.t J = 2^^ 



(b) Error w.r.t e 



Figure 6: errp£(t = 0.13) for AP scheme 
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Figure 8: errpe(t = 0.05) for time splitting scheme 
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Figure 9: errje(t = 0.05) for time splitting scheme 




Figure 10: eiTpe{t = 0.13) for time splitting scheme 
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Figure 12: Re('u^) at T = 0.05 for e = 0.005 
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(a) J = 128 



(b) J = 256 



Figure 13: Re(ti^) at T = 0.13 for e = 0.005 




(a) J = 512 

Figure 14: Re('u'^) at T = 0.13 for e = 0.005 



3.3. Numerical experiments in dimension 2. The previous subsection was de- 
voted to one dimensional simulations. We present here the equivalent error anal- 
ysis for the two dimensional case. The computational domain is now the square 
[—0.5, 1.5]^ discretized with J points in each direction. We still need a reference 
solution generated with Strang splitting scheme applied to equation (1.1)-(1.2). We 
take Jref = = 8192, so Ax ~ 2.5 10"^ and At = e/100. In a first test, the initial 
datum is related to the one chosen for the one dimensional case 

-25 



vo{x) 



1 



Vln(e^" + e-5n 



where r = 
values of J 



(X — 0.5)^ + (y — 0.5)^. We evaluate the error function for different 
2*^, where M is chosen in [4, 11], and various scaled Planck parameter 
e in [5.10~^, 10*^]. The formation time of singularities is as in one dimension situation 
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located around 0.10. The reference particle and current densities for e = 5.10 



are represented on Figures 15 and 17 before and after formation of singularities. 
We clearly see that a tiny front appears after shocks both in particle and current 
densities. The wave function presents strong oscillations with new ones created past 



the shock creation (see Fig. 17) 



o.g| 

0.8 



- 0.5 I 
0-4 I 
0.3 I 
0.2 I 



o 



I 



' 0, 



(a) pr, 




Figure 15: Contour plot of particle and current densities at time t = 0.05 and for 
e = 5.10-"^ 




The error analysis results are equivalent to the ones obtained for one dimensional 
simulation. We recover an independence of the error with respect to e before the 
formation of singularities (see Figs 18 and 19). We always need to take finer mesh 
past the shocks formation (see Figs 20 and 21). 
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The last computation concerns a non-symmetric solution. The aim of our simu- 
lations is to show that we recover good qualitative properties, even for non radially 
symmetric data. We consider here as initial datum a Maxwell distribution with two 
temperatures 9i = 0.05 and 62 = 0.015 for the initial amplitude given by 

1 1 / (x- 0.5)2 (y-0.5f 

«o = —p^ / / exp 



The initial phase is reduced to 0o = (note that from (1.11 ), dt4i\t=o — "I'^op 7^ 0, 
so (j) becomes instantaneously non-trivial). The scaled Planck factor e is equal 
to e = 0.005 and the simulations are performed with J = 2048 intervals in each 
direction. To distinguish the evolution for each direction {x,y), we make contour 
plots of wave function and physical observables but also present slice for plane 
X = 0.5 and y = 0.5 respectively on bottom right and upper right axes. The 
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Figure 20: err pe(t = 0.13) for AP scheme 



solution has many oscillations in both directions before formation of singularities 
(see the real part of the wave function on Fig. 22 at time t = 0.035). The 
oscillatory nature is enhanced at time t = 0.08 mainly in y— direction (see Fig. 23). 



We recover the good behavior of physical observables before the formation time 
of singularities. Actually, we compare the particle densities computed by Strang 
splitting scheme and our AP scheme at time t = 0.035 on Figures 24 and 25 The 
similar representation for is available on Figures 28 and [29j Past the formation 
time of singularities, the contour plots could be thought as equivalent both for 
particle and current densities. The only noticeable differences can be seen on x = 0.5 
and y = 0.5 plane shces (see Fig. [26) |27| |30] and |3l]) . 
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Figure 23: Re{u^){t = 0.08) 
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Figure 26: Strang = 0.08) 
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Figure 29: AP = 0.035) 




Figure 31: AP = 0.08) 
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4. Extension to other frameworks 



4.1. The hnear case: global solution of the viscous eikonal equation. The 

advantage of introducing an artificial diffusion is more striking in the linear case. 
Consider 

^2 



iedtu'^ H Au^ 

2 



i4>o/e 



(4.1) 



with T4xt = K;xt(ii x) real- valued. As noticed in |22j, using the Madelung transform 
in this context is interesting only if vacuum can be avoided. Having the idea of 
Grenier [29j in mind, the counterpart of (1.8) reads 



i-Aa', 



cxt 



0, 



*|i=0 



We note that the two equations are decoupled: the second equation is of Hamilton- 
Jacobi type (known as the eikonal equation in this context), and the solution must 
be expected to develop singularities in finite time (see e.g. [E]), so this approach is 
doomed to fail for large time. On the other hand, introducing the artificial diffusion 
as in (1.11), we obtain 



z-Aa^ 
2 



1, 



dt<j>' + ^\Vcl)'\-' + V, 



ext 



''|i=0 



t=0 



(4.2) 



The system is still decoupled, but the good news is that the diffusion introduced into 
the Hamilton- Jacobi equation smoothes the solution out. Therefore, this system can 
be considered as a good candidate to obtain an AP scheme for the linear Schrodinger 
equation, regardless of the presence of vacuum. 



Instead of giving full details of the analogue of Theorems 1.4 1.6 and 1.7 as well 
as of their proofs, we shall simply give a functional framework where the viscous 
eikonal equation can by solved globally in time. Since we consider the case e > 
only, we drop out the e"^ term in the viscous eikonal equation, and consider 



For ^ 1, let 

SLfc = {/ G C^(M'^ 

For ^ 2, let 

SQfc = {/ G C\M!' 
The key result is the following: 



+ Kxt = A(/)cik; (/'eik(0, x) = (l)o{x). 

;M), 5"/ G L°°(M'^), VI ^ |aK A:}. 
;M), 5"/ G L°°(M'^), V2 < |a| ^ A:}. 



(4.3) 



Lemma 4.1. Let k ^ 2 and (po G SLk, V^ext 
unique solution 0cik G C7(M+; SLfc^i). 



^loc 



(M+;SLfc). Then (4.3) has a 



Remark 4.2. No such global result must be expected in the larger class SQfc. Indeed 
suppose d = 1, Vext = and (^(^{x) — ""'^ 
C(M+; QL2), then w = d'^(peik e C(M+; L 

dtw + vdxW + up' = d'^w; w{0, x) = —1. 



-X /2. If (4.3) had a solution 0cik £ 
would solve 
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The solution of this equation does not depend on x. It is given by w{t, x) = 1/ (t—l), 
hence a contradiction. However, in the periodic case x G T'^, this technical issue 
disappears (see Remark |4.3| below) . 



Proof. The first step consists in constructing the gradient of (j)eik- Differentiating 



(4.3) in space, we have to solve 

dtv + v-Vv + VVc^t = ^v; v{0, x) = \/(j)o{x). (4.4) 

This equation is solved locally in time by a fixed point argument by using Duhamel's 
formula 

v{t) = e*^V0o - r e^'-'^^{v ■ Vv){s)ds - f e^'-'^^{VVe^t){s)ds. 



Jo Jo 
The right hand side is a contraction on 

{weC ([0,T]; VF^-i'°°) , ||w;||Lco([o,T];W'fc-i.-) ^ 2||V0o|Ih/^-i.-} , 

provided that T > is sufficiently small. This follows from the properties of the 
heat kernel: for t > 0, 

l|e*^/IU- ^ II/IIl-; l|e*^v/|Uoo ^ ^II/IIl-. 



The solution to (4.4) is then global, thanks to the maximum principle (see e.g. |47| 
Proposition 52.8], which implies that u_ ^ v ^ where 

dtv± =F ||Vycxt||Ls° = 0; ^±(0,2;) = ±||V(/)o||l-'- 
Once such a solution v G C(M+; W'^~^'°°) is constructed, set 

4>cik{t) = 4>o- j Q|w(s)P + Vext{s) - divf(s)^ ds. 
We check that 5t(V0eik — v) = Vdt(peik — dtv = 0. Therefore, v = V(/>eik, and 



solves (4.3). Rewriting (4.3) as 



dt(t>eik + ^I^^P + V^xt = A^eik; </'eik(0, x) = 4>o{x), 

Duhamel's formula reads 

0eik(t) = e*^0o - fe'^'-'^'^{\v{s)\')ds- f e^'-'^''V,,t{s)ds, 
Jo Jo 

from which we conclude that (/>eik G C(M-(_; SLfc_i). □ 

Remark 4.3. In the periodic setting x G T*^ = (R/(27rZ))'^, it is natural to work in 
the Wiener algebra 

W=lf:T''^C, fix) = J2 h^'''" ^ith {bj)^^^, G £^(Z« 

and for k ^ 0, 

W'' = {f ■.T'^^ C, a"/ G W, V|a| ^ A;}. 



Then Lemma 4.1 is easily adapted by replacing SL^ with W^. 
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4.2. Other nonlinearities. One might beheve that Theorem |1.7| is bomid to the 
cubic one-dimensional Schrodinger equation, which is completely integrable, and 



that the two aspects are related. We show that this is not the case, since Theorem 1.7 
remains valid for other nonlinearities (with d = 1). Consider now 

=.2 

iedtu' + jAu' = f u^ (i, G M+ x R'^. (4.5) 

We suppose that there exists 6 > such that f'{y) ^ 6, for all y ^ 0. 

Following [29], the only modification we have to make to recover the results in 



Section 2.1 consists in replacing the symmetrizer (|2.2|) with 







From our assumption on /, S and its inverse are uniformly bounded, provided 
that a is bounded. As a matter of fact, the exact assumption in [29j is /' > 0, and 



all the results in Section 2.1 remain valid under this assumption. 

The more precise assumption f'{y) ^ 5 becomes useful to prove that in dimension 
d = 1, the solution is global. The main point to notice is that the conservation of 



mass and momentum are the same as for (1.1), like in Proposition |1.1[ On the other 



hand, the conservation of energy becomes 

I (ikV^^(t)lli. + £ F {\u%t, dx^ = 0, 
where F is an anti-derivative of /, 



F{y) = I f{r)dr. 

By assumption. 



y 







F{y) ^ + f{o)y, 

so the potential energy controls the L^-norm: 

h%t)\\l, [ F{\u%t,x)\')dx-fm\u%t)\\l, 



6 
2 

< - 
6 



F{\u%t,x)\')dx-fmu%0)\\ 



where we have used the conservation of mass in the last inequality. This implies an 
estimate of the form 

\\eVu^t)\\l, + \\n^t)\\l^C, 
with C independent of t ^ and e G (0,1]. Therefore, the analysis presented in 



Section 2.3 can be repeated line by line. 



Example 4.4. The assumption /' ^ (5 > is satisfied in the following cases: 

• Cubic-quintic nonlinearity: f{y) = y + \y^, A ^ 0. 

• Cubic plus saturated nonlinearity: f{y) = 6y + V i^xy ^ ^^Vi^ > 0- 
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